Fit VBGE models to back-calculated length-at-age

Load libraries

pkgs <- c("tidyverse", "tidylog", "rTPC", "nls.multstart", "broom", "RColorBrewer", "ggridges",
          "viridis", "minpack.lm", "patchwork", "ggtext", "brms", "tidybayes", "modelr", "investr")

# minpack.lm needed if using nlsLM()
if(length(setdiff(pkgs, rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  
  }

invisible(lapply(pkgs, library,character.only = T))

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.2     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Attaching package: 'tidylog'


The following objects are masked from 'package:dplyr':

    add_count, add_tally, anti_join, count, distinct, distinct_all,
    distinct_at, distinct_if, filter, filter_all, filter_at, filter_if,
    full_join, group_by, group_by_all, group_by_at, group_by_if,
    inner_join, left_join, mutate, mutate_all, mutate_at, mutate_if,
    relocate, rename, rename_all, rename_at, rename_if, rename_with,
    right_join, sample_frac, sample_n, select, select_all, select_at,
    select_if, semi_join, slice, slice_head, slice_max, slice_min,
    slice_sample, slice_tail, summarise, summarise_all, summarise_at,
    summarise_if, summarize, summarize_all, summarize_at, summarize_if,
    tally, top_frac, top_n, transmute, transmute_all, transmute_at,
    transmute_if, ungroup


The following objects are masked from 'package:tidyr':

    drop_na, fill, gather, pivot_longer, pivot_wider, replace_na,
    spread, uncount


The following object is masked from 'package:stats':

    filter


Loading required package: viridisLite

Loading required package: Rcpp

Loading 'brms' package (version 2.19.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').


Attaching package: 'brms'


The following object is masked from 'package:stats':

    ar



Attaching package: 'tidybayes'


The following objects are masked from 'package:brms':

    dstudent_t, pstudent_t, qstudent_t, rstudent_t


The following objects are masked from 'package:ggridges':

    scale_point_color_continuous, scale_point_color_discrete,
    scale_point_colour_continuous, scale_point_colour_discrete,
    scale_point_fill_continuous, scale_point_fill_discrete,
    scale_point_size_continuous



Attaching package: 'modelr'


The following object is masked from 'package:broom':

    bootstrap

# devtools::install_github("seananderson/ggsidekick") ## not on CRAN 
library(ggsidekick)
theme_set(theme_sleek())

# Load functions
home <- here::here()
fxn <- list.files(paste0(home, "/R/functions"))
invisible(sapply(FUN = source, paste0(home, "/R/functions/", fxn)))

Load cache

# qwraps2::lazyload_cache_dir(path = paste0(home, "/R/analyze-data/02-fit-vbge_cache/html"))

Read and trim data

d <- #read_csv(paste0(home, "/data/for-analysis/dat.csv")) |> 
  read_csv("https://raw.githubusercontent.com/maxlindmark/perch-growth/master/data/for-analysis/dat.csv") |> 
  filter(age_ring == "Y", # use only length-at-age by filtering on age_ring
         !area %in% c("VN", "TH")) 

Rows: 338460 Columns: 12
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (6): age_ring, area, gear, ID, sex, keep
dbl (6): length_mm, age_bc, age_catch, catch_year, cohort, final_length

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
filter: removed 43,886 rows (13%), 294,574 rows remaining

# Sample size by area and cohort
ns <- d |> 
  group_by(cohort, area) |> 
  summarise(n = n())

group_by: 2 grouping variables (cohort, area)
summarise: now 431 rows and 3 columns, one group variable remaining (cohort)

# Minimum number of observations per area and cohort
d$area_cohort <- as.factor(paste(d$area, d$cohort))

d <- d |>
  group_by(area_cohort) |> 
  filter(n() > 100)

group_by: one grouping variable (area_cohort)
filter (grouped): removed 2,637 rows (1%), 291,937 rows remaining

# Minimum number of observations per area, cohort, and age
d$area_cohort_age <- as.factor(paste(d$area, d$cohort, d$age_bc))

d <- d |>
  group_by(area_cohort_age) |> 
  filter(n() > 10)

group_by: one grouping variable (area_cohort_age)
filter (grouped): removed 3,521 rows (1%), 288,416 rows remaining

# Minimum number of cohorts in a given area
cnt <- d |>
  group_by(area) |>
  summarise(n = n_distinct(cohort)) |>
  filter(n >= 10)

group_by: one grouping variable (area)
summarise: now 10 rows and 2 columns, ungrouped
filter: no rows removed

d <- d[d$area %in% cnt$area, ]

# Plot cleaned data
ggplot(d, aes(age_bc, length_mm)) +
  geom_jitter(size = 0.1, height = 0, alpha = 0.1) +
  scale_x_continuous(breaks = seq(20)) +
  theme(axis.text.x = element_text(angle = 0)) +
  theme(axis.text = element_text(size = 12), axis.title = element_text(size = 15)) +
  labs(x = "Age", y = "Length (mm)") +
  guides(color = "none") + 
  facet_wrap(~area, scale = "free_x")

# Longitude and latitude attributes for each area
area <- c("BS", "BT", "FB", "FM", "HO", "JM", "MU", "RA", "SI_EK", "SI_HA", "TH")
nareas <- length(area)
lat <- c(60, 60.4, 60.3, 60.5, 63.7, 58, 59, 65.9, 57.3, 57.4, 56.1)
lon <- c(21.5, 18.1, 19.5, 18, 20.9, 16.8, 18.1, 22.3, 16.6, 16.7, 15.9)
area_attr <- data.frame(cbind(area = area, lat = lat, lon = lon)) |>
  mutate_at(c("lat", "lon"), as.numeric)

mutate_at: converted 'lat' from character to double (0 new NA)
           converted 'lon' from character to double (0 new NA)

Fit VBGE models

# Get individual growth parameters (functions: VBGF/Gompertz,nls_out,fit_nls)
IVBG <- d |>
  group_by(ID) |>
  summarize(nls_out(fit_nls(length_mm, age_bc, min_nage = 5, model = "VBGF")))

group_by: one grouping variable (ID)

summarize: now 75,245 rows and 5 columns, ungrouped

Inspect predictions

IVBG <- IVBG |> drop_na(k) # The NAs are because the model didn't converge or because they were below the threshold age

drop_na: removed 51,640 rows (69%), 23,605 rows remaining

IVBG <- IVBG |>
  mutate(k_lwr = k - 1.96*k_se,
         k_upr = k + 1.96*k_se,
         linf_lwr = linf - 1.96*linf_se,
         linf_upr = linf + 1.96*linf_se,
         row_id = row_number())

mutate: new variable 'k_lwr' (double) with 23,603 unique values and 0% NA

        new variable 'k_upr' (double) with 23,603 unique values and 0% NA

        new variable 'linf_lwr' (double) with 23,603 unique values and 0% NA

        new variable 'linf_upr' (double) with 23,603 unique values and 0% NA

        new variable 'row_id' (integer) with 23,605 unique values and 0% NA

# Plot all K's
IVBG |>
  #filter(row_id < 5000) |>
  ggplot(aes(row_id, k, ymin = k_lwr, ymax = k_upr)) +
  geom_point(alpha = 0.2) +
  geom_errorbar(alpha = 0.2) +
  NULL

# Plot all L_inf's
IVBG |>
  #filter(row_id < 5000) |>
  ggplot(aes(row_id, linf, ymin = linf_lwr, ymax = linf_upr)) +
  geom_point(alpha = 0.2) +
  geom_errorbar(alpha = 0.2) +
  NULL

# We can also consider removing individuals where the SE of k is larger than the fit
IVBG |>
  #mutate(keep = ifelse(k > quantile(k_se, probs = 0.95), "N", "Y")) |>
  mutate(keep = ifelse(k < k_se, "N", "Y")) |>
  #filter(row_id < 10000) |>
  ggplot(aes(row_id, k, ymin = k_lwr, ymax = k_upr, color = keep)) +
  geom_point(alpha = 0.5) +
  facet_wrap(~keep, ncol = 1) +
  geom_errorbar(alpha = 0.5) +
  NULL

mutate: new variable 'keep' (character) with 2 unique values and 0% NA

# Add back cohort and area variables
IVBG <- IVBG |> 
  left_join(d |> select(ID, area_cohort)) |> 
  separate(area_cohort, into = c("area", "cohort"), remove = TRUE, sep = " ") |> 
  mutate(cohort = as.numeric(cohort))

Adding missing grouping variables: `area_cohort_age`

select: dropped 11 variables (length_mm, age_bc, age_catch, age_ring, area, …)

Joining with `by = join_by(ID)`
left_join: added 2 columns (area_cohort_age, area_cohort)
> rows only in x 0
> rows only in y (146,393)
> matched rows 142,023 (includes duplicates)
> =========
> rows total 142,023
mutate: converted 'cohort' from character to double (0 new NA)

# Summarise and save for sample size
sample_size <- IVBG |>
  group_by(area) |> 
  summarise(n_cohort = length(unique(cohort)),
            min_cohort = min(cohort),
            max_cohort = min(cohort),
            n_individuals = length(unique(ID)),
            n_data_points = n())

group_by: one grouping variable (area)
summarise: now 10 rows and 6 columns, ungrouped

sample_size

# A tibble: 10 × 6
   area  n_cohort min_cohort max_cohort n_individuals n_data_points
   <chr>    <int>      <dbl>      <dbl>         <int>         <int>
 1 BS          17       1985       1985          1334          7688
 2 BT          22       1972       1972           961          5371
 3 FB          47       1969       1969          6040         37381
 4 FM          37       1963       1963          5057         32642
 5 HO          29       1982       1982          1150          6210
 6 JM          60       1953       1953          4868         28749
 7 MU          18       1981       1981          1104          6666
 8 RA          14       1985       1985           572          3122
 9 SI_EK       24       1968       1968           658          3571
10 SI_HA       38       1967       1967          1861         10623

sample_size |>
  ungroup() |>
  summarise(sum_ind = sum(n_individuals), sum_n = sum(n_data_points))

ungroup: no grouping variables
summarise: now one row and 2 columns, ungrouped

# A tibble: 1 × 2
  sum_ind  sum_n
    <int>  <int>
1   23605 142023

write_csv(sample_size, paste0(home, "/output/sample_size.csv"))

# Compare how the means and quantiles differ depending on this filtering
IVBG_filter <- IVBG |>
  drop_na(k_se) |>
  #filter(k_se < quantile(k_se, probs = 0.95)) |> 
  filter(k_se < k)

drop_na: no rows removed
filter: removed 4,192 rows (3%), 137,831 rows remaining

# Summarize growth coefficients by cohort and area
VBG <- IVBG |>
  filter(k_se < k) |> # new!
  group_by(cohort, area) |>
  summarize(k_mean = mean(k, na.rm = T),
            k_median = quantile(k, prob = 0.5, na.rm = T),
            linf_median = quantile(linf, prob = 0.5, na.rm = T),
            k_lwr = quantile(k, prob = 0.05, na.rm = T),
            k_upr = quantile(k, prob = 0.95, na.rm = T)) |> 
  ungroup()

filter: removed 4,192 rows (3%), 137,831 rows remaining
group_by: 2 grouping variables (cohort, area)
summarize: now 306 rows and 7 columns, one group variable remaining (cohort)
ungroup: no grouping variables

VBG_filter <- IVBG_filter |>
  group_by(cohort, area) |>
  summarize(k_mean = mean(k, na.rm = T),
            k_median = quantile(k, prob = 0.5, na.rm = T),
            k_lwr = quantile(k, prob = 0.05, na.rm = T),
            k_upr = quantile(k, prob = 0.95, na.rm = T)) |> 
  ungroup()

group_by: 2 grouping variables (cohort, area)
summarize: now 306 rows and 6 columns, one group variable remaining (cohort)
ungroup: no grouping variables

ggplot() +
  geom_ribbon(data = VBG, aes(cohort, k_median, ymin = k_lwr, ymax = k_upr,
                              fill = "All k's"), alpha = 0.4) +
  geom_ribbon(data = VBG_filter, aes(cohort, k_median, ymin = k_lwr, ymax = k_upr,
                                     fill = "Filtered"), alpha = 0.4) +
  geom_line(data = VBG, aes(cohort, k_median, color = "All k's")) + 
  geom_line(data = VBG_filter, aes(cohort, k_median, color = "Filtered"), linetype = 2) + 
  guides(fill = "none") +
  facet_wrap(~area)

ggplot() +
  geom_line(data = VBG, aes(cohort, k_median, color = "All k's")) + 
  geom_line(data = VBG_filter, aes(cohort, k_median, color = "Filtered"), linetype = 2) + 
  guides(fill = "none") +
  facet_wrap(~area)

# No difference at all really. We should probably just discuss that with this model, achieving biologically reasonable parameter values and a good fit to data are sometimes two different things. In our case, we just want a representative value of the growth (as in time to reach average maximum size in the population) that accounts for the entire growth trajectory of an individual, for each area and cohort.

BT 1996 has a very high k, so I’ll inspect it in more detail

Add GAM-predicted temperature to growth models

pred_temp <- read_csv(paste0(home, "/output/gam_predicted_temps.csv")) |> 
  rename(cohort = year)

Rows: 380 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): area, model
dbl (2): year, temp

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed one variable (cohort)

VBG <- VBG |>
  left_join(pred_temp, by = c("area", "cohort"))

left_join: added 2 columns (temp, model)
           > rows only in x     0
           > rows only in y  ( 74)
           > matched rows     306
           >                 =====
           > rows total       306

# Save data for map-plot
cohort_sample_size <- IVBG |>
  group_by(area, cohort) |> 
  summarise(n = n()) # individuals, not samples!

group_by: 2 grouping variables (area, cohort)
summarise: now 306 rows and 3 columns, one group variable remaining (area)

VBG <- left_join(VBG, cohort_sample_size, by = c("cohort", "area"))

left_join: added one column (n)
           > rows only in x     0
           > rows only in y  (  0)
           > matched rows     306
           >                 =====
           > rows total       306

write_csv(VBG, paste0(home, "/output/vbg.csv"))

# Calculate the plotting order (also used for map plot)
order <- VBG |>
  ungroup() |>
  group_by(area) |>
  summarise(min_temp = min(temp)) |>
  arrange(desc(min_temp))

ungroup: no grouping variables
group_by: one grouping variable (area)
summarise: now 10 rows and 2 columns, ungrouped

order

# A tibble: 10 × 2
   area  min_temp
   <chr>    <dbl>
 1 SI_HA     7.83
 2 BS        6.22
 3 BT        5.87
 4 FM        5.87
 5 SI_EK     5.49
 6 MU        5.16
 7 FB        5.04
 8 HO        3.99
 9 JM        3.37
10 RA        3.06

write_csv(order, paste0(home, "/output/ranked_temps.csv"))

nareas <- length(unique(order$area)) + 2 # to skip the brightest colors that are hard to read
colors <- colorRampPalette(brewer.pal(name = "RdYlBu", n = 10))(nareas)[-c(6,7)]

Plot VBGE fits

# Sample 30 IDs per area and plot their data and VBGE fits
set.seed(4)
ids <- IVBG |> distinct(ID, .keep_all = TRUE) |> group_by(area) |> slice_sample(n = 30)

distinct: removed 118,418 rows (83%), 23,605 rows remaining

group_by: one grouping variable (area)

slice_sample (grouped): removed 23,305 rows (99%), 300 rows remaining

#ids |> ungroup() |> group_by(area) |> summarise(n = length(unique(ID))) |> arrange(n)

IVBG2 <- IVBG |>
  filter(ID %in% ids$ID) |> 
  distinct(ID, .keep_all = TRUE) |> 
  select(ID, k, linf)

filter: removed 140,319 rows (99%), 1,704 rows remaining

distinct: removed 1,404 rows (82%), 300 rows remaining

select: dropped 10 variables (k_se, linf_se, k_lwr, k_upr, linf_lwr, …)

d2 <- d |>
  ungroup() |>
  filter(ID %in% ids$ID) |>
  left_join(IVBG2, by = "ID") |>
  mutate(length_mm_pred = linf*(1-exp(-k*age_bc)))

ungroup: no grouping variables

filter: removed 286,712 rows (99%), 1,704 rows remaining

left_join: added 2 columns (k, linf)

           > rows only in x       0

           > rows only in y  (    0)

           > matched rows     1,704

           >                 =======

           > rows total       1,704

mutate: new variable 'length_mm_pred' (double) with 1,697 unique values and 0% NA

ggplot(d2, aes(age_bc, length_mm, group = ID, color = ID)) +
  geom_jitter(width = 0.3, height = 0, alpha = 0.5, size = 0.4) +
  geom_line(aes(age_bc, length_mm_pred, group = ID, color = ID),
            inherit.aes = FALSE, alpha = 0.8, linewidth = 0.3) + 
  guides(color = "none") +
  scale_color_viridis(discrete = TRUE, name = "Area", option = "cividis") +
  labs(x = "Age (years)", y = "Length (mm)") +
  facet_wrap(~factor(area, levels = (arrange(area_attr, desc(lat)))$area), ncol = 4)

ggsave(paste0(home, "/figures/vb_fits.pdf" ), width = 17, height = 12, units = "cm")

k <- IVBG |> 
  ggplot(aes(factor(area, order$area), k, fill = factor(area, order$area))) + 
  coord_cartesian(ylim = c(0, 0.7)) +
  geom_violin(alpha = 0.8, color = NA) +  
  geom_boxplot(outlier.color = NA, width = 0.2, alpha = 0.9, fill = NA, size = 0.4) +
  scale_fill_manual(values = colors, name = "Area") +
  guides(fill = "none", color = "none") +
  labs(x = "Area", y = expression(italic(k)))

k

linf <- IVBG |> 
  filter(linf < 2300) |> 
  filter(linf > 130) |> 
  ggplot(aes(factor(area, order$area), linf, fill = factor(area, order$area))) + 
  geom_violin(alpha = 0.8, color = NA) +  
  geom_boxplot(outlier.color = NA, width = 0.1, alpha = 0.9, fill = NA, size = 0.4) +
  coord_cartesian(ylim = c(0, 2000)) +
  scale_fill_manual(values = colors, name = "Area") +
  guides(fill = "none", color = "none") +
  labs(x = "", y = expression(paste(italic(L[infinity]), " [mm]")))

filter: removed 1,218 rows (1%), 140,805 rows remaining

filter: no rows removed

linf

k / linf

ggsave(paste0(home, "/figures/vb_pars.pdf" ), width = 17, height = 17, units = "cm")

Fit Sharpe-Schoolfield model to k

By area

model <- 'sharpeschoolhigh_1981'

# Get starting values on full dataset for Sharpe-Schoolfield model
dat <- VBG |>
  select(k_median, temp, area) |>
  rename(rate = k_median)

select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)

rename: renamed one variable (rate)

lower <- get_lower_lims(dat$temp, dat$rate, model_name = model)
upper <- get_upper_lims(dat$temp, dat$rate, model_name = model)
start <- get_start_vals(dat$temp, dat$rate, model_name = model)
  
# Sharpe-Schoolfield model fit to data for each area
preds <- NULL
pred <- NULL
pars <- list()
t_opts <- list()

for(a in unique(dat$area)) {
  
  # Get data
  dd <- dat |> filter(area == a)
  
  # Fit model
  fit <- nls_multstart(
    rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 8),
    data = dd,
    iter = c(3, 3, 3, 3),
    start_lower = start*0.5,
    start_upper = start*2,
    lower = lower,
    upper = upper,
    supp_errors = 'Y'
    )
  
  # Make predictions on new data
  new_data <- data.frame(temp = seq(min(dd$temp), max(dd$temp), length.out = 100))
  
  pred <- augment(fit, newdata = new_data) |> mutate(area = a)
  
  # Add 95% CI's
  
  if(unique(dd$area) %in% c("MU", "HO")){
    
    pred$lwr <- NA
    pred$upr <- NA
    
  } else{
  
    pred$upr <- as_tibble(predFit(fit, newdata = new_data, interval = "confidence", level = 0.95))$upr
    pred$lwr <- as_tibble(predFit(fit, newdata = new_data, interval = "confidence", level = 0.95))$lwr  
    
  }
    
  # Add to general data frame
  preds <- data.frame(rbind(preds, pred))
  
  # Add parameter estimates
  pars[[a]] <- summary(fit)$coefficients |>
    as.data.frame() |>
    rownames_to_column(var = "parameter") |>
    mutate(area = a)
  
  # Extract t_topt
  t_opts[[a]] <- data.frame(t_opt = get_topt(fit),
                            area = a)

}

filter: removed 246 rows (80%), 60 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 269 rows (88%), 37 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 268 rows (88%), 38 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 282 rows (92%), 24 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 259 rows (85%), 47 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 284 rows (93%), 22 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 288 rows (94%), 18 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 277 rows (91%), 29 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 289 rows (94%), 17 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 292 rows (95%), 14 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

# For comparison with the bayesian model below
nls_pars <- bind_rows(pars)
t_opts <- bind_rows(t_opts)

nls_pars <- nls_pars |> 
  mutate(upr = Estimate + 1.96 *`Std. Error`,
         lwr = Estimate - 1.96 *`Std. Error`) |> 
  dplyr::select(parameter, Estimate, upr, lwr, area) |> 
  mutate(source = "nls")

mutate: new variable 'upr' (double) with 40 unique values and 0% NA

        new variable 'lwr' (double) with 40 unique values and 0% NA

mutate: new variable 'source' (character) with one unique value and 0% NA

All areas pooled

# One sharpe schoolfield fitted to all data
fit_all <- nls_multstart(
    rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 8),
    data = dat,
    iter = c(3, 3, 3, 3),
    start_lower = start*0.5,
    start_upper = start*2,
    lower = lower,
    upper = upper,
    supp_errors = 'Y'
    )

nls_pooled_pars <- summary(fit_all)$coefficients |>
  as.data.frame() |>
  rownames_to_column(var = "parameter") |> 
  dplyr::select(parameter, Estimate) |> 
  pivot_wider(names_from = parameter, values_from = Estimate)

pivot_wider: reorganized (parameter, Estimate) into (r_tref, e, eh, th) [was 4x2, now 1x4]

# Make predictions on new data
new_data_all <- data.frame(temp = seq(min(dat$temp), max(dat$temp), length.out = 100))

pred_all <- augment(fit_all, newdata = new_data_all) |> 
  mutate(area = "all")

mutate: new variable 'area' (character) with one unique value and 0% NA

# To get confidence intervals around our nls fits, we can do use "investr"
# https://stackoverflow.com/questions/61341287/how-to-calculate-confidence-intervals-for-nonlinear-least-squares-in-r
# https://stackoverflow.com/questions/52973626/how-to-calculate-95-prediction-interval-from-nls

pred_all$upr <- as_tibble(predFit(fit_all, newdata = new_data_all, interval = "confidence", level = 0.95))$upr
pred_all$lwr <- as_tibble(predFit(fit_all, newdata = new_data_all, interval = "confidence", level = 0.95))$lwr
  
# Add t_opt (not correct equation I think!) from Padfield 2021 ISME
kb <- 8.62e-05
# data.frame(par = names(coef(fit_all)), est = coef(fit_all)) |> 
#   pivot_wider(names_from = par, values_from = est) |> 
#   summarise(t_opt = (eh*th) / (eh + kb*th*log((eh/e)-1)))
# 
# get_topt

# This rTPC function just finds the temperature where the function is maximized, I use the same approach brms fits
get_topt(fit_all)

[1] 9.741584

# Bootstrapping uncertainty... ? Need to refit I think
fit_nlsLM <- minpack.lm::nlsLM(rate~sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 8),
                               data = dat,
                               start = start,
                               lower = start*0.5,
                               upper = start*2)

# bootstrap using case resampling
library(car)

Loading required package: carData

Attaching package: 'car'

The following object is masked from 'package:dplyr':

    recode

The following object is masked from 'package:purrr':

    some

boot1 <- Boot(fit_nlsLM, method = 'case', R = 10000)

Loading required namespace: boot

head(boot1$t)

        r_tref         e        eh       th
[1,] 0.3645396 0.2819862 0.6132375 9.641969
[2,] 0.3689331 0.2819862 0.6132375 9.094924
[3,] 0.3755589 0.2819862 0.6132375 9.377010
[4,] 0.3599203 0.2819862 0.6132375 9.731394
[5,] 0.3808942 0.2819862 0.6132375 9.374603
[6,] 0.3759139 0.2819862 0.6132375 9.152401

str(head(boot1$t))

 num [1:6, 1:4] 0.365 0.369 0.376 0.36 0.381 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:4] "r_tref" "e" "eh" "th"

boot1$t |> 
  as_tibble() |> 
  pivot_longer(everything()) |> 
  ggplot(aes(value)) + 
  geom_histogram() + 
  facet_wrap(~name, scales = "free")

pivot_longer: reorganized (r_tref, e, eh, th) into (name, value) [was 10000x4, now 40000x2]

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Calculate t_opt by adding predictions
# https://padpadpadpad.github.io/rTPC/articles/bootstrapping_models.html
boot1_preds <- boot1$t |>
  as.data.frame() |>
  drop_na() |>
  mutate(iter = 1:n()) |>
  group_by_all() |>
  do(data.frame(temp = seq(min(dat$temp), max(dat$temp), length.out = 1000))) |>
  ungroup() |>
  mutate(pred = sharpeschoolhigh_1981(temp, r_tref, e, eh, th, tref = 8))

drop_na: no rows removed
mutate: new variable 'iter' (integer) with 10,000 unique values and 0% NA
group_by_all: 5 grouping variables (r_tref, e, eh, th, iter)
ungroup: no grouping variables
mutate: new variable 'pred' (double) with 10,000,000 unique values and 0% NA

# Find the maximum for each boostrapped parameter combination ("iter")
boot1_preds |> 
  group_by(iter) |> 
  filter(pred == max(pred)) |> 
  ggplot(aes(temp)) + 
  geom_histogram()

group_by: one grouping variable (iter)
filter (grouped): removed 9,990,000 rows (>99%), 10,000 rows remaining
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Plot Sharpe Schoolfield fits

# Plot growth coefficients by year and area against mean SST

p1 <- preds |>
  ggplot(aes(temp, .fitted, color = factor(area, order$area))) + 
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 1, alpha = 0.8) +
  geom_line(linewidth = 1) +
  geom_line(data = pred_all, aes(temp, .fitted), linewidth = 1, inherit.aes = FALSE, linetype = 2) +
  scale_color_manual(values = colors, name = "Area") +
  guides(color = guide_legend(nrow = 1, reverse = TRUE, title.position = "top", title.hjust = 0.5,
                              override.aes = list(size = 1))) +
  scale_x_continuous(breaks = seq(-5, 20, 2)) +
  labs(x = "Temperature (°C)",
       y = "von Bertalanffy growth coefficient (*k*)") +
  theme(axis.title.y = ggtext::element_markdown(),
        legend.position = "bottom",
        legend.direction = "horizontal")

# Now facet and add ribbons
p1 +
  facet_wrap(~factor(area, levels = (arrange(area_attr, desc(lat)))$area)) + 
  geom_ribbon(aes(ymin = lwr, ymax = upr, fill = factor(area, order$area)), alpha = 0.3, color = NA) +
  geom_ribbon(data = pred_all, aes(temp, .fitted, ymin = lwr, ymax = upr), alpha = 0.3, color = NA) +
  scale_fill_manual(values = colors, name = "Area") +
  guides(fill = "none") +
  geom_line(linewidth = 1)

Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
-Inf

Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
-Inf

ggsave(paste0(home, "/figures/supp/sharpe_school_nls.pdf" ), width = 17, height = 17, units = "cm")

Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
-Inf

Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
-Inf

Can we fit a single Sharpe Schoolfield with area-specific parameters with brms?

# Again, here are the data we are fitting:
ggplot(dat, aes(temp, rate, color = factor(area, order$area))) + 
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 0.6)  

# Which family to use? Well we can try with a Gaussian or a student-t first... 
# After that, maybe gamma or lognormal, if we predict negative K's
hist(dat$rate)

# Here's the equation:
sharpeschoolhigh_1981

function (temp, r_tref, e, eh, th, tref) 
{
    tref <- 273.15 + tref
    k <- 8.62e-05
    boltzmann.term <- r_tref * exp(e/k * (1/tref - 1/(temp + 
        273.15)))
    inactivation.term <- 1/(1 + exp(eh/k * (1/(th + 273.15) - 
        1/(temp + 273.15))))
    return(boltzmann.term * inactivation.term)
}
<bytecode: 0x7ff102716828>
<environment: namespace:rTPC>

# Add in fixed parameters
dat$bk <- 8.62e-05
dat$tref <- 8 + 273.15

# We can use the nls() parameters as starting values
avg_nls_pars <- nls_pars |> 
  group_by(parameter) |> 
  summarise(mean_par = mean(Estimate)) |> 
  pivot_wider(names_from = parameter, values_from = mean_par)

group_by: one grouping variable (parameter)

summarise: now 4 rows and 2 columns, ungrouped

pivot_wider: reorganized (parameter, mean_par) into (e, eh, r_tref, th) [was 4x2, now 1x4]

avg_nls_pars

# A tibble: 1 × 4
      e    eh r_tref    th
  <dbl> <dbl>  <dbl> <dbl>
1 0.611  2.75  0.309  13.0

# As a first guess, I use the nls estimate as the mean, and standard deviation something close to it. We can also bound som fo them by 0, since negative values shouldn't be possible
# (better visualization including bounds further down)
n=10000
hist(rnorm(avg_nls_pars$e, 0.5, n=n), main = "e", xlim = c(0, 5))

hist(rnorm(avg_nls_pars$eh, 1, n=n), main = "eh", xlim = c(0, 5))

hist(rnorm(avg_nls_pars$r_tref, 0.5, n=n), main = "r_tref", xlim = c(0, 2.5))

hist(rnorm(avg_nls_pars$th, 5, n=n), main = "th", xlim = c(0, 30))

# Set these priors
prior <- c(prior(normal(0.6112299, 0.5), nlpar = "e", lb = 0),
           prior(normal(2.753837, 1), nlpar = "eh", lb = 0),
           prior(normal(0.3085462, 0.5), nlpar = "rtref", lb = 0),
           prior(normal(13.03523, 5), nlpar = "th", lb = 0))

fit_prior <- brm(
  bf(rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))),
     rtref + e + eh + th ~ 1, 
     nl = TRUE),
  data = dat,
  cores = 2,
  chains = 2,
  iter = 1500,
  sample_prior = "only",
  seed = 9,
  prior = prior
)

Compiling Stan program...

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

Start sampling

Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

# Global prior predictive check in relation to data. Doesn't loo too informative...
dat |>
  data_grid(temp = seq_range(temp, n = 51)) |> 
  ungroup() |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) |> 
  add_epred_draws(fit_prior) |> 
  mutate(r_tref = nls_pooled_pars$r_tref, e = nls_pooled_pars$e, eh = nls_pooled_pars$eh, th = nls_pooled_pars$th) |> 
  mutate(man_pred = (r_tref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15))))) |> 
  ggplot(aes(temp, y = .epred)) +
  stat_lineribbon(.width = c(0.95), alpha = 0.3, color = "black", fill = "black") + 
  geom_line(aes(y = man_pred), color = "tomato3", linetype = 2) + 
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 0.6) +
  scale_color_manual(values = colors, name = "Area") +
  NULL

ungroup: no grouping variables

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

mutate (grouped): new variable 'r_tref' (double) with one unique value and 0% NA

                  new variable 'e' (double) with one unique value and 0% NA

                  new variable 'eh' (double) with one unique value and 0% NA

                  new variable 'th' (double) with one unique value and 0% NA

mutate (grouped): new variable 'man_pred' (double) with 51 unique values and 0% NA

# Now make sure it converges with real data but no random effects
fit_global <- brm(
  bf(rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))),
     rtref + e + eh + th ~ 1,
     nl = TRUE),
  data = dat,
  cores = 4,
  chains = 4,
  iter = 4000, 
  seed = 9,
  sample_prior = "yes",
  prior = prior,
  control = list(adapt_delta = 0.99, max_treedepth = 12)
)

Compiling Stan program...

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

Start sampling

plot(fit_global)

pp_check(fit_global)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

plot(conditional_effects(fit_global, effect = "temp"), points = TRUE)

dat |>
  data_grid(temp = seq_range(temp, n = 51)) |> 
  ungroup() |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) |> 
  add_epred_draws(fit_global) |> 
  # 
  mutate(r_tref = nls_pooled_pars$r_tref, e = nls_pooled_pars$e, eh = nls_pooled_pars$eh, th = nls_pooled_pars$th) |>
  mutate(man_pred = (r_tref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15))))) |> 
  ggplot(aes(temp, y = .epred)) +
  stat_lineribbon(.width = c(0.95), alpha = 0.3, color = "black", fill = "black") + 
  geom_line(aes(y = man_pred), color = "tomato3", linetype = 2) + 
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 0.6) +
  scale_color_manual(values = colors, name = "Area")

ungroup: no grouping variables

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

mutate (grouped): new variable 'r_tref' (double) with one unique value and 0% NA

                  new variable 'e' (double) with one unique value and 0% NA

                  new variable 'eh' (double) with one unique value and 0% NA

                  new variable 'th' (double) with one unique value and 0% NA

mutate (grouped): new variable 'man_pred' (double) with 51 unique values and 0% NA

# Plot prior vs posterior
post <- fit_global |>
  as_draws_df() |>
  dplyr::select(b_rtref_Intercept, b_e_Intercept, b_eh_Intercept, b_th_Intercept) |> 
  rename(r_tref = b_rtref_Intercept, 
         e = b_e_Intercept,
         eh = b_eh_Intercept,
         th = b_th_Intercept) |> 
  pivot_longer(everything(), names_to = "parameter") |> 
  mutate(type = "Posterior")

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (r_tref, e, eh, th)

pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

prior_draws <- fit_global |>
  as_draws_df() |>
  dplyr::select(prior_b_rtref, prior_b_e, prior_b_eh, prior_b_th) |> 
  rename(r_tref = prior_b_rtref, 
         e = prior_b_e,
         eh = prior_b_eh,
         th = prior_b_th) |> 
  pivot_longer(everything(), names_to = "parameter") |> 
  mutate(type = "Prior")

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (r_tref, e, eh, th)

pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

dist <- bind_rows(prior_draws, post)

ggplot(dist, aes(value, fill = type)) +
  geom_density(color = NA, alpha = 0.5) +
  scale_fill_brewer(palette = "Set1") +
  facet_wrap(~parameter, scales = "free") + 
  theme(legend.position = c(0.35, 0.9)) +
  labs(fill = "")

Ok, seems to work with priors and everything. They are roughly as broad as can be and still have a fit. Now fit the full model, with random area effects, more iterations and chains!

# Still good enough, might change some priors. Now look at random area effects!
# Gaussian model
fitb <- brm(
  bf(rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))),
     rtref + e + eh + th ~ 1 + (1|area),
     nl = TRUE),
  data = dat,
  cores = 4,
  chains = 4,
  iter = 5000,
  sample_prior = "yes",
  save_pars = save_pars(all = TRUE),
  seed = 9,
  prior = prior,
  control = list(adapt_delta = 0.99, max_treedepth = 12)
)

Compiling Stan program...

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

Start sampling

# Student-t model
fitbs <- brm(
  bf(rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))),
     rtref + e + eh + th ~ 1 + (1|area),
     nl = TRUE),
  data = dat,
  family = student,
  cores = 4,
  chains = 4,
  iter = 5000,
  sample_prior = "yes",
  save_pars = save_pars(all = TRUE),
  seed = 9,
  prior = prior,
  control = list(adapt_delta = 0.99, max_treedepth = 12)
)

Compiling Stan program...
Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

Start sampling

print(fitbs, digits = 4)

 Family: student 
  Links: mu = identity; sigma = identity; nu = identity 
Formula: rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15))))/(1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))) 
         rtref ~ 1 + (1 | area)
         e ~ 1 + (1 | area)
         eh ~ 1 + (1 | area)
         th ~ 1 + (1 | area)
   Data: dat (Number of observations: 306) 
  Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
         total post-warmup draws = 10000

Group-Level Effects: 
~area (Number of levels: 10) 
                    Estimate Est.Error l-95% CI u-95% CI   Rhat Bulk_ESS
sd(rtref_Intercept)   0.0322    0.0232   0.0020   0.0899 1.0007     1732
sd(e_Intercept)       0.1879    0.1283   0.0109   0.4904 1.0008     2638
sd(eh_Intercept)      0.4462    0.4463   0.0150   1.5712 1.0006     2023
sd(th_Intercept)      1.1362    1.0236   0.0454   3.6889 1.0002     3409
                    Tail_ESS
sd(rtref_Intercept)     2690
sd(e_Intercept)         3115
sd(eh_Intercept)        4198
sd(th_Intercept)        4052

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI   Rhat Bulk_ESS Tail_ESS
rtref_Intercept   0.3152    0.1168   0.1962   0.6317 1.0012     1404     2871
e_Intercept       0.5867    0.3309   0.0879   1.3781 1.0010     1627     2523
eh_Intercept      1.7546    0.6191   0.7373   3.3783 1.0025     2292     1961
th_Intercept     12.2022    4.1609   4.7838  20.7330 1.0011     1529     2959

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI   Rhat Bulk_ESS Tail_ESS
sigma   0.0397    0.0021   0.0357   0.0438 1.0009     9550     7205
nu     24.7640   13.4838   7.9340  59.1615 1.0002     9928     8129

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Compare fitted models based on ELPD. Preferred model in the first row!

fitb_loo <- loo(fitb, moment_match = TRUE)

Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.

fitbs_loo <- loo(fitbs, moment_match = TRUE)

loo_compare(fitb_loo, fitbs_loo)

      elpd_diff se_diff
fitbs  0.0       0.0   
fitb  -0.2       1.3   

Quick test to see if we can use a different set of parameters, e.g., “flat” (normal but wide and bounded) but tight

n=10000
hist(rnorm(avg_nls_pars$e, 1, n=n), main = "e", xlim = c(0, 2))

hist(rnorm(avg_nls_pars$eh, 3, n=n), main = "eh", xlim = c(0, 5))

hist(rnorm(avg_nls_pars$r_tref, 0.5, n=n), main = "r_tref", xlim = c(0, 1))

hist(rnorm(avg_nls_pars$th, 15, n=n), main = "th", xlim = c(3, 20))

prior_flat <- c(prior(normal(0.6112299, 1), nlpar = "e", lb = 0, ub = 2),
                prior(normal(2.753837, 3), nlpar = "eh", lb = 0, ub = 5),
                prior(normal(0.3085462, 0.5), nlpar = "rtref", lb = 0, ub = 1),
                prior(normal(13.03523, 15), nlpar = "th", lb = 3, ub = 20))

# Student-t model with "flat" priors
fitbs_flat <- brm(
  bf(rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))),
     rtref + e + eh + th ~ 1 + (1|area),
     nl = TRUE),
  data = dat,
  family = student,
  cores = 4,
  chains = 4,
  iter = 6000,
  sample_prior = "yes",
  save_pars = save_pars(all = TRUE),
  seed = 9,
  prior = prior_flat,
  control = list(adapt_delta = 0.99, max_treedepth = 12)
)

Compiling Stan program...

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

Start sampling

print(fitbs_flat, digits = 4)

 Family: student 
  Links: mu = identity; sigma = identity; nu = identity 
Formula: rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15))))/(1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))) 
         rtref ~ 1 + (1 | area)
         e ~ 1 + (1 | area)
         eh ~ 1 + (1 | area)
         th ~ 1 + (1 | area)
   Data: dat (Number of observations: 306) 
  Draws: 4 chains, each with iter = 6000; warmup = 3000; thin = 1;
         total post-warmup draws = 12000

Group-Level Effects: 
~area (Number of levels: 10) 
                    Estimate Est.Error l-95% CI u-95% CI   Rhat Bulk_ESS
sd(rtref_Intercept)   0.0395    0.0313   0.0020   0.1182 1.0005     2315
sd(e_Intercept)       0.1880    0.1311   0.0096   0.5056 1.0006     3308
sd(eh_Intercept)      0.3715    0.3886   0.0141   1.3542 1.0031     2393
sd(th_Intercept)      1.0269    0.9518   0.0315   3.4494 1.0012     3386
                    Tail_ESS
sd(rtref_Intercept)     4256
sd(e_Intercept)         4662
sd(eh_Intercept)        4541
sd(th_Intercept)        4968

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI   Rhat Bulk_ESS Tail_ESS
rtref_Intercept   0.3849    0.1621   0.1994   0.8170 1.0035     1642     2433
e_Intercept       0.6985    0.4325   0.0862   1.7128 1.0031     1730     4029
eh_Intercept      1.5619    0.6696   0.3278   3.2472 1.0010     3031     2471
th_Intercept     10.4526    4.4138   3.7388  19.1941 1.0039     1622     3638

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI   Rhat Bulk_ESS Tail_ESS
sigma   0.0397    0.0020   0.0358   0.0437 1.0001    11322     8949
nu     24.7114   13.5015   8.0481  59.5055 1.0001    12506    10073

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

print(fitbs, digits = 4)

 Family: student 
  Links: mu = identity; sigma = identity; nu = identity 
Formula: rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15))))/(1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))) 
         rtref ~ 1 + (1 | area)
         e ~ 1 + (1 | area)
         eh ~ 1 + (1 | area)
         th ~ 1 + (1 | area)
   Data: dat (Number of observations: 306) 
  Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
         total post-warmup draws = 10000

Group-Level Effects: 
~area (Number of levels: 10) 
                    Estimate Est.Error l-95% CI u-95% CI   Rhat Bulk_ESS
sd(rtref_Intercept)   0.0322    0.0232   0.0020   0.0899 1.0007     1732
sd(e_Intercept)       0.1879    0.1283   0.0109   0.4904 1.0008     2638
sd(eh_Intercept)      0.4462    0.4463   0.0150   1.5712 1.0006     2023
sd(th_Intercept)      1.1362    1.0236   0.0454   3.6889 1.0002     3409
                    Tail_ESS
sd(rtref_Intercept)     2690
sd(e_Intercept)         3115
sd(eh_Intercept)        4198
sd(th_Intercept)        4052

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI   Rhat Bulk_ESS Tail_ESS
rtref_Intercept   0.3152    0.1168   0.1962   0.6317 1.0012     1404     2871
e_Intercept       0.5867    0.3309   0.0879   1.3781 1.0010     1627     2523
eh_Intercept      1.7546    0.6191   0.7373   3.3783 1.0025     2292     1961
th_Intercept     12.2022    4.1609   4.7838  20.7330 1.0011     1529     2959

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI   Rhat Bulk_ESS Tail_ESS
sigma   0.0397    0.0021   0.0357   0.0438 1.0009     9550     7205
nu     24.7640   13.4838   7.9340  59.1615 1.0002     9928     8129

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

The student model is preferred. Inspect it!

# Check fit
summary(fitbs)

 Family: student 
  Links: mu = identity; sigma = identity; nu = identity 
Formula: rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15))))/(1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))) 
         rtref ~ 1 + (1 | area)
         e ~ 1 + (1 | area)
         eh ~ 1 + (1 | area)
         th ~ 1 + (1 | area)
   Data: dat (Number of observations: 306) 
  Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
         total post-warmup draws = 10000

Group-Level Effects: 
~area (Number of levels: 10) 
                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(rtref_Intercept)     0.03      0.02     0.00     0.09 1.00     1732     2690
sd(e_Intercept)         0.19      0.13     0.01     0.49 1.00     2638     3115
sd(eh_Intercept)        0.45      0.45     0.02     1.57 1.00     2023     4198
sd(th_Intercept)        1.14      1.02     0.05     3.69 1.00     3409     4052

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rtref_Intercept     0.32      0.12     0.20     0.63 1.00     1404     2871
e_Intercept         0.59      0.33     0.09     1.38 1.00     1627     2523
eh_Intercept        1.75      0.62     0.74     3.38 1.00     2292     1961
th_Intercept       12.20      4.16     4.78    20.73 1.00     1529     2959

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.04      0.00     0.04     0.04 1.00     9550     7205
nu       24.76     13.48     7.93    59.16 1.00     9928     8129

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(fitbs)

pp_check(fitbs) + theme(legend.position = c(0.1, 0.9))

Using 10 posterior draws for ppc type 'dens_overlay' by default.

ggsave(paste0(home, "/figures/supp/sharpe_brms_ppcheck.pdf" ), width = 11, height = 11, units = "cm")

# Plot prior vs posterior
post <- fitbs |>
  as_draws_df() |>
  dplyr::select(b_rtref_Intercept, b_e_Intercept, b_eh_Intercept, b_th_Intercept) |> 
  rename(r_tref = b_rtref_Intercept, 
         e = b_e_Intercept,
         eh = b_eh_Intercept,
         th = b_th_Intercept) |> 
  pivot_longer(everything(), names_to = "parameter") |> 
  mutate(type = "Posterior")

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (r_tref, e, eh, th)

pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 10000x4, now 40000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

prior_draws <- fitbs |>
  as_draws_df() |>
  dplyr::select(prior_b_rtref, prior_b_e, prior_b_eh, prior_b_th) |> 
    rename(r_tref = prior_b_rtref, 
         e = prior_b_e,
         eh = prior_b_eh,
         th = prior_b_th) |> 
  pivot_longer(everything(), names_to = "parameter") |> 
  mutate(type = "Prior")

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (r_tref, e, eh, th)

pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 10000x4, now 40000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

dist <- bind_rows(prior_draws, post)

ggplot(dist, aes(value, fill = type)) +
  geom_density(color = NA, alpha = 0.5) +
  scale_fill_brewer(palette = "Set1") +
  facet_wrap(~parameter, scales = "free") + 
  theme(legend.position = c(0.35, 0.9)) +
  labs(fill = "")

ggsave(paste0(home, "/figures/supp/sharpe_brms_prior.pdf" ), width = 17, height = 17, units = "cm")

# Now do the "flat" prior model
# Plot prior vs posterior
post <- fitbs_flat |>
  as_draws_df() |>
  dplyr::select(b_rtref_Intercept, b_e_Intercept, b_eh_Intercept, b_th_Intercept) |> 
  rename(r_tref = b_rtref_Intercept, 
         e = b_e_Intercept,
         eh = b_eh_Intercept,
         th = b_th_Intercept) |> 
  pivot_longer(everything(), names_to = "parameter") |> 
  mutate(type = "Posterior")

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (r_tref, e, eh, th)

pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 12000x4, now 48000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

prior_draws <- fitbs_flat |>
  as_draws_df() |>
  dplyr::select(prior_b_rtref, prior_b_e, prior_b_eh, prior_b_th) |> 
    rename(r_tref = prior_b_rtref, 
         e = prior_b_e,
         eh = prior_b_eh,
         th = prior_b_th) |> 
  pivot_longer(everything(), names_to = "parameter") |> 
  mutate(type = "Prior")

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (r_tref, e, eh, th)

pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 12000x4, now 48000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

dist <- bind_rows(prior_draws, post)

ggplot(dist, aes(value, fill = type)) +
  geom_density(color = NA, alpha = 0.5) +
  scale_fill_brewer(palette = "Set1") +
  facet_wrap(~parameter, scales = "free") + 
  theme(legend.position = c(0.35, 0.9)) +
  labs(fill = "")

Calculate T_opt

ts <- dat |> 
  data_grid(temp = seq_range(temp, n = 100)) |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) |> 
  add_epred_draws(fitbs, re_formula = NA, ndraws = 10000) 

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

ggplot(ts) + 
  geom_line(aes(temp, y = .epred, group = .draw), alpha = .1)

# Compute quantiles across spaghetties
t_opt <- ts |> 
  group_by(.draw) |> 
  filter(.epred == max(.epred)) |> 
  ungroup() |> 
  filter(!temp == min(temp)) |> # remove values where there was no optimum
  filter(!temp == max(temp))

group_by: one grouping variable (.draw)

filter (grouped): removed 990,000 rows (99%), 10,000 rows remaining

ungroup: no grouping variables

filter: removed 87 rows (1%), 9,913 rows remaining

filter: removed 213 rows (2%), 9,700 rows remaining

quantile(t_opt$temp, probs = c(0.025, 0.5, 0.975))

     2.5%       50%     97.5% 
 6.402944  9.468607 14.206451 

ggplot(t_opt, aes(temp)) + 
  geom_histogram(bins = length(unique(t_opt$temp))) + 
  theme_void()

# To add the histogram inside the main plot, I can convert it to "numbers-at-temperature" and plot as a bin. This way I know the maximum X. I will plot these as separatly. In ggplot, you cannot technically do this, but you can scale the axes to have similar values, the change labels.
t_opt2 <- t_opt |> 
  group_by(temp) |> 
  summarise(n = n())

group_by: one grouping variable (temp)

summarise: now 98 rows and 2 columns, ungrouped

ggplot(t_opt2, aes(temp, n)) + 
  geom_bar(stat = "identity", color = NA, fill = "grey30", width = 0.14)

Warning: `position_stack()` requires non-overlapping x intervals

scaleFactor <- max(t_opt2$n) / max(dat$rate)

# Make the main plot (conditional effect of temperature, with and without random effects)
# Predictions without random effects
glob <- dat |>
  data_grid(temp = seq_range(temp, n = 100)) |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) |> 
  add_epred_draws(fitbs, re_formula = NA) |> 
  ungroup()

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

ungroup: no grouping variables

dat |>
  group_by(area) |>
  data_grid(temp = seq_range(temp, n = 100)) |> 
  ungroup() |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) |> 
  add_epred_draws(fitbs) |> 
  ungroup() |> 
  ggplot(aes(temp, y = .epred, color = factor(area, order$area))) +
  stat_lineribbon(data = glob, aes(temp, .epred), .width = c(0.9), inherit.aes = FALSE,
                 fill = "black", color = "black", alpha = 0.1) +
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 1, alpha = 0.8) +
  stat_lineribbon(.width = c(0)) +
  stat_lineribbon(data = glob, aes(temp, .epred), .width = c(0), inherit.aes = FALSE,
                 color = "black", alpha = 0.9, linetype = 2) +
  #coord_cartesian(ylim = c(min(dat$rate) - 0.01, sort(dat$rate)[length(dat$rate) - 1] + 0.01)) +
  coord_cartesian(expand = 0) +
  scale_color_manual(values = colors, name = "Area") +
  guides(fill = "none",
         color = guide_legend(nrow = 1, reverse = TRUE, title.position = "top", title.hjust = 0.5,
                              override.aes = list(size = 1))) +
  theme(axis.title.y = ggtext::element_markdown(),
        legend.position = "bottom",
        legend.direction = "horizontal") +
  labs(x = "Temperature (°C)",
       y = "von Bertalanffy growth coefficient (*k*)") + 
  geom_bar(data = t_opt2, aes(temp, n/scaleFactor/4), stat = "identity", inherit.aes = FALSE,
           color = NA, fill = "gray70", alpha = 0.9, width = 0.16) + 
  scale_y_continuous(sec.axis = sec_axis(~.*scaleFactor*4, name = "", breaks = NULL, labels = NULL)) +
  NULL

group_by: one grouping variable (area)
ungroup: no grouping variables

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

ungroup: no grouping variables

Warning: `position_stack()` requires non-overlapping x intervals

ggsave(paste0(home, "/figures/sharpe_school_bayes.pdf" ), width = 17, height = 17, units = "cm")

Warning: `position_stack()` requires non-overlapping x intervals

Area-specific T_opts as a ridgeplot

# Calculate T_opt by area
tsa <- dat |> 
  data_grid(area = unique(dat$area),
            temp = seq_range(temp, n = 100)) |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) |> 
  add_epred_draws(fitbs, re_formula = NULL, ndraws = 10000) 

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

t_opt_a <- tsa |> 
  group_by(.draw, area) |> 
  filter(.epred == max(.epred)) |> 
  ungroup() |> 
  filter(!temp == min(temp)) |> # remove values where there was no optimum
  filter(!temp == max(temp))

group_by: 2 grouping variables (.draw, area)

filter (grouped): removed 9,900,000 rows (99%), 100,000 rows remaining

ungroup: no grouping variables

filter: removed 6,652 rows (7%), 93,348 rows remaining

filter: removed 8,255 rows (9%), 85,093 rows remaining

t_opt_a_sum <- t_opt_a |> 
  group_by(area) |> 
  summarise(median = median(temp),
            lwr = quantile(temp, probs = 0.1),
            upr = quantile(temp, probs = 0.9))

group_by: one grouping variable (area)

summarise: now 10 rows and 4 columns, ungrouped

ggplot(t_opt, aes(temp)) + 
  geom_histogram(bins = length(unique(t_opt$temp)))

# Compute quantiles across spaghetties
rect_dat <- data.frame(xmin = c(quantile(t_opt$temp, probs = c(0.05)), quantile(t_opt$temp, probs = c(0.1))),
                       xmax = c(quantile(t_opt$temp, probs = c(0.95)), quantile(t_opt$temp, probs = c(0.9))),
                       ymin = c(-Inf, -Inf),
                       ymax = c(Inf, Inf),
                       interval = c("90%", "80%"))

ggplot(t_opt_a, aes(temp, factor(area, levels = order$area), fill = factor(area, levels = order$area))) + 
  geom_vline(xintercept = quantile(t_opt$temp, probs = c(0.5)), alpha = 0.6, linetype = 2, linewidth = 0.75) +
  geom_rect(data = rect_dat |> filter(interval == "90%"), aes(xmin = xmin, xmax = xmax), ymin = -Inf, ymax = Inf, 
            inherit.aes = FALSE, fill = "grey", alpha = 0.5) +
  geom_rect(data = rect_dat |> filter(interval == "80%"), aes(xmin = xmin, xmax = xmax), ymin = -Inf, ymax = Inf, 
            inherit.aes = FALSE, fill = "grey", alpha = 0.5) +
  geom_density_ridges(alpha = 0.7, color = NA, scale = 1,
                      stat = "binline", bins = 50) +
  scale_fill_manual(values = colors, name = "Area") +
  scale_color_manual(values = colors, name = "Area") +
  labs(x = "Temperature (°C)", y = "") +
  guides(fill = "none") +
  geom_errorbar(data = t_opt_a_sum, aes(xmin = lwr, xmax = upr, y = area, color = factor(area, levels = order$area)), linewidth = 1,
                inherit.aes = FALSE, width = 0) +
  geom_point(data = t_opt_a_sum, aes(median, area, fill = factor(area, levels = order$area)), shape = 21, color = "white", size = 3) +
  NULL

filter: removed one row (50%), one row remaining

filter: removed one row (50%), one row remaining

ggsave(paste0(home, "/figures/t_opt_ridges.pdf" ), width = 17, height = 17, units = "cm")


# Add the nls t_opts
ggplot(t_opt_a, aes(temp, factor(area, levels = order$area), fill = factor(area, levels = order$area))) + 
  geom_vline(xintercept = quantile(t_opt$temp, probs = c(0.5)), alpha = 0.6, linetype = 2, linewidth = 0.75) +
  geom_rect(data = rect_dat |> filter(interval == "90%"), aes(xmin = xmin, xmax = xmax), ymin = -Inf, ymax = Inf, 
            inherit.aes = FALSE, fill = "grey", alpha = 0.5) +
  geom_rect(data = rect_dat |> filter(interval == "80%"), aes(xmin = xmin, xmax = xmax), ymin = -Inf, ymax = Inf, 
            inherit.aes = FALSE, fill = "grey", alpha = 0.5) +
  geom_density_ridges(alpha = 0.8, color = NA, scale = 0.9,
                      stat = "binline", bins = 50) +
  scale_fill_manual(values = colors, name = "Area") +
  #scale_color_manual(values = colors, name = "Area") +
  labs(x = "Temperature (°C)", y = "") +
  guides(fill = "none") +
  geom_point(data = t_opts, aes(t_opt, area, fill = area), shape = 21, color = "white", size = 5) +
  NULL

filter: removed one row (50%), one row remaining
filter: removed one row (50%), one row remaining

Plot comparisons with the nls multstart fits

# Compare with area-specific sharpe scool!
area_pred_brms <- dat |>
  group_by(area) |> 
  data_grid(temp = seq_range(temp, n = 100)) |> 
  ungroup() |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) |> 
  add_epred_draws(fitbs)

group_by: one grouping variable (area)

ungroup: no grouping variables

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

# Plot area specific in comparison with nls multstart
p1 <- ggplot(area_pred_brms, aes(temp, y = .epred, color = factor(area, order$area),
                           fill = factor(area, order$area))) +
  stat_lineribbon(.width = c(0.95), alpha = 0.4, color = NA) +
  stat_lineribbon(.width = c(0)) + 
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 0.6) +
  scale_color_manual(values = colors, name = "Area") +
  scale_fill_manual(values = colors, name = "Area") +
  scale_linetype_manual(values = 2) +
  facet_wrap(~factor(area, rev(order$area))) +
  guides(fill = "none",
         color = guide_legend(nrow = 1, reverse = TRUE, title.position = "top", title.hjust = 0.5,
                              override.aes = list(size = 1))) +
  theme(axis.title.y = ggtext::element_markdown(),
        legend.position = "bottom",
        legend.direction = "horizontal") +
  labs(x = "Temperature (°C)",
       y = "von Bertalanffy growth coefficient (*k*)",
       linetype = "")

p1

ggsave(paste0(home, "/figures/supp/sharpe_school_bayes_ci_facet.pdf" ), width = 17, height = 17, units = "cm")

# Now compare to nls fits
p1 + geom_line(data = preds, aes(temp, .fitted, linetype = "nls multstart"), color = "gray10")

# And for the population-level prediction
nd_brms <- dat |>
  data_grid(temp = seq_range(temp, n = 100)) |> 
  ungroup() |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15)

ungroup: no grouping variables

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

pred_brms <- predict(fitbs, newdata = nd_brms, re_formula = NA) |> as.data.frame()
nd_brms$pred_pop <- pred_brms$Estimate

# (why not also comparing the global brms model...)
pred_brms_global <- predict(fit_global, newdata = nd_brms) |> as.data.frame()
nd_brms$pred_glob <- pred_brms_global$Estimate

ggplot(nd_brms, aes(temp, y = pred_pop, color = "mixed brms model")) +
  geom_point(data = dat, aes(temp, rate), size = 0.6, color = "grey30") +
  geom_line() + 
  geom_line(data = pred_all, aes(temp, .fitted, color = "nls"), linewidth = 1, inherit.aes = FALSE) +
  geom_line(data = nd_brms, aes(temp, pred_glob, color = "global brms model"), linewidth = 1, inherit.aes = FALSE) +
  guides(fill = "none",
         color = guide_legend(ncol = 1, reverse = TRUE, title.position = "top", title.hjust = 0.5,
                              override.aes = list(size = 1)),
         linetype = guide_legend(keywidth = unit(1, "cm"))) +
  theme(axis.title.y = ggtext::element_markdown()) +
  labs(x = "Temperature (°C)",
       y = "von Bertalanffy growth coefficient (*k*)",
       color = "")

And a comparison of parameter estimates…

get_variables(fitbs)

  [1] "b_rtref_Intercept"              "b_e_Intercept"                 
  [3] "b_eh_Intercept"                 "b_th_Intercept"                
  [5] "sd_area__rtref_Intercept"       "sd_area__e_Intercept"          
  [7] "sd_area__eh_Intercept"          "sd_area__th_Intercept"         
  [9] "sigma"                          "nu"                            
 [11] "r_area__rtref[BS,Intercept]"    "r_area__rtref[BT,Intercept]"   
 [13] "r_area__rtref[FB,Intercept]"    "r_area__rtref[FM,Intercept]"   
 [15] "r_area__rtref[HO,Intercept]"    "r_area__rtref[JM,Intercept]"   
 [17] "r_area__rtref[MU,Intercept]"    "r_area__rtref[RA,Intercept]"   
 [19] "r_area__rtref[SI_EK,Intercept]" "r_area__rtref[SI_HA,Intercept]"
 [21] "r_area__e[BS,Intercept]"        "r_area__e[BT,Intercept]"       
 [23] "r_area__e[FB,Intercept]"        "r_area__e[FM,Intercept]"       
 [25] "r_area__e[HO,Intercept]"        "r_area__e[JM,Intercept]"       
 [27] "r_area__e[MU,Intercept]"        "r_area__e[RA,Intercept]"       
 [29] "r_area__e[SI_EK,Intercept]"     "r_area__e[SI_HA,Intercept]"    
 [31] "r_area__eh[BS,Intercept]"       "r_area__eh[BT,Intercept]"      
 [33] "r_area__eh[FB,Intercept]"       "r_area__eh[FM,Intercept]"      
 [35] "r_area__eh[HO,Intercept]"       "r_area__eh[JM,Intercept]"      
 [37] "r_area__eh[MU,Intercept]"       "r_area__eh[RA,Intercept]"      
 [39] "r_area__eh[SI_EK,Intercept]"    "r_area__eh[SI_HA,Intercept]"   
 [41] "r_area__th[BS,Intercept]"       "r_area__th[BT,Intercept]"      
 [43] "r_area__th[FB,Intercept]"       "r_area__th[FM,Intercept]"      
 [45] "r_area__th[HO,Intercept]"       "r_area__th[JM,Intercept]"      
 [47] "r_area__th[MU,Intercept]"       "r_area__th[RA,Intercept]"      
 [49] "r_area__th[SI_EK,Intercept]"    "r_area__th[SI_HA,Intercept]"   
 [51] "prior_b_rtref"                  "prior_b_e"                     
 [53] "prior_b_eh"                     "prior_b_th"                    
 [55] "prior_sigma"                    "prior_nu"                      
 [57] "prior_sd_area"                  "prior_sd_area__1"              
 [59] "prior_sd_area__2"               "prior_sd_area__3"              
 [61] "lprior"                         "lp__"                          
 [63] "z_1[1,1]"                       "z_1[1,2]"                      
 [65] "z_1[1,3]"                       "z_1[1,4]"                      
 [67] "z_1[1,5]"                       "z_1[1,6]"                      
 [69] "z_1[1,7]"                       "z_1[1,8]"                      
 [71] "z_1[1,9]"                       "z_1[1,10]"                     
 [73] "z_2[1,1]"                       "z_2[1,2]"                      
 [75] "z_2[1,3]"                       "z_2[1,4]"                      
 [77] "z_2[1,5]"                       "z_2[1,6]"                      
 [79] "z_2[1,7]"                       "z_2[1,8]"                      
 [81] "z_2[1,9]"                       "z_2[1,10]"                     
 [83] "z_3[1,1]"                       "z_3[1,2]"                      
 [85] "z_3[1,3]"                       "z_3[1,4]"                      
 [87] "z_3[1,5]"                       "z_3[1,6]"                      
 [89] "z_3[1,7]"                       "z_3[1,8]"                      
 [91] "z_3[1,9]"                       "z_3[1,10]"                     
 [93] "z_4[1,1]"                       "z_4[1,2]"                      
 [95] "z_4[1,3]"                       "z_4[1,4]"                      
 [97] "z_4[1,5]"                       "z_4[1,6]"                      
 [99] "z_4[1,7]"                       "z_4[1,8]"                      
[101] "z_4[1,9]"                       "z_4[1,10]"                     
[103] "accept_stat__"                  "stepsize__"                    
[105] "treedepth__"                    "n_leapfrog__"                  
[107] "divergent__"                    "energy__"                      

summary(fitbs)

 Family: student 
  Links: mu = identity; sigma = identity; nu = identity 
Formula: rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15))))/(1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))) 
         rtref ~ 1 + (1 | area)
         e ~ 1 + (1 | area)
         eh ~ 1 + (1 | area)
         th ~ 1 + (1 | area)
   Data: dat (Number of observations: 306) 
  Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
         total post-warmup draws = 10000

Group-Level Effects: 
~area (Number of levels: 10) 
                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(rtref_Intercept)     0.03      0.02     0.00     0.09 1.00     1732     2690
sd(e_Intercept)         0.19      0.13     0.01     0.49 1.00     2638     3115
sd(eh_Intercept)        0.45      0.45     0.02     1.57 1.00     2023     4198
sd(th_Intercept)        1.14      1.02     0.05     3.69 1.00     3409     4052

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rtref_Intercept     0.32      0.12     0.20     0.63 1.00     1404     2871
e_Intercept         0.59      0.33     0.09     1.38 1.00     1627     2523
eh_Intercept        1.75      0.62     0.74     3.38 1.00     2292     1961
th_Intercept       12.20      4.16     4.78    20.73 1.00     1529     2959

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.04      0.00     0.04     0.04 1.00     9550     7205
nu       24.76     13.48     7.93    59.16 1.00     9928     8129

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

pop_pars <- summary(fitbs)$fixed |> 
  as.data.frame() |> 
  rownames_to_column(var = "parameter") |> 
  mutate(parameter = str_remove(parameter, "_Intercept"),
         parameter = ifelse(parameter == "rtref", "r_tref", parameter)) |> 
  rename(pop_Estimate = Estimate) |> 
  dplyr::select(parameter, pop_Estimate)

mutate: changed 4 values (100%) of 'parameter' (0 new NA)

rename: renamed one variable (pop_Estimate)

pop_pars <- fitbs |> 
  gather_draws(b_rtref_Intercept,
               b_e_Intercept,
               b_eh_Intercept,
               b_th_Intercept) |> 
  rename(parameter = .variable) |> 
  mutate(parameter = str_remove(parameter, "b_"),
         parameter = str_remove(parameter, "_Intercept"),
         parameter = ifelse(parameter == "rtref", "r_tref", parameter)) |> 
  rename(.pop_value = .value )

rename: renamed one variable (parameter)

mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)

rename: renamed one variable (.pop_value)

brms_pars <- fitbs |> 
  gather_draws(r_area__rtref[area, Intercept],
               r_area__e[area, Intercept],
               r_area__eh[area, Intercept],
               r_area__th[area, Intercept]) |> 
  rename(parameter = .variable) |> 
  mutate(parameter = str_remove(parameter, "r_area__"),
         parameter = ifelse(parameter == "rtref", "r_tref", parameter)) |> 
  dplyr::select(-Intercept) |>
  left_join(pop_pars, by = c(".chain", ".iteration", ".draw", "parameter")) |> 
  mutate(.value = .pop_value + .value) |> 
  summarise(Estimate = median(.value),
            lwr = quantile(.value, probs = 0.025),
            upr = quantile(.value, probs = 0.975)) |> 
  ungroup() |> 
  mutate(source = "brms")

rename: renamed one variable (parameter)

mutate (grouped): changed 400,000 values (100%) of 'parameter' (0 new NA)

Adding missing grouping variables: `Intercept`
left_join: added one column (.pop_value)
> rows only in x 0
> rows only in y ( 0)
> matched rows 400,000
> =========
> rows total 400,000
mutate (grouped): changed 400,000 values (100%) of '.value' (0 new NA)
summarise: now 40 rows and 6 columns, 2 group variables remaining (area,
Intercept)
ungroup: no grouping variables
mutate: new variable 'source' (character) with one unique value and 0% NA

pars <- bind_rows(brms_pars, nls_pars)

# Need to trim errorbars...
ggplot(pars, aes(area, Estimate, color = source)) + 
  geom_point() +
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0) + 
  facet_wrap(~parameter, scales = "free")

pars |> 
  ggplot(aes(area, Estimate, color = source)) + 
  geom_point(position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0,
                position = position_dodge(width = 0.5)) + 
  coord_cartesian(ylim = c(-3, 20)) + 
  facet_wrap(~parameter, scales = "free")

# And without errorbars...
pars |> 
  ggplot(aes(area, Estimate, color = source)) + 
  geom_point(position = position_dodge(width = 0.5)) +
  facet_wrap(~parameter, scales = "free")

pars |> as.data.frame()

    area Intercept parameter     Estimate           lwr          upr source
1     BS Intercept         e  0.601959304  6.223733e-02 1.491306e+00   brms
2     BS Intercept        eh  1.562468135  3.569021e-01 3.681207e+00   brms
3     BS Intercept    r_tref  0.272467401  1.769527e-01 6.117746e-01   brms
4     BS Intercept        th 11.695436560  4.398969e+00 2.110572e+01   brms
5     BT Intercept         e  0.605091810  7.097435e-02 1.504592e+00   brms
6     BT Intercept        eh  1.594603464  6.644355e-01 3.669490e+00   brms
7     BT Intercept    r_tref  0.294454470  1.997231e-01 6.628758e-01   brms
8     BT Intercept        th 12.521423712  4.883776e+00 2.150019e+01   brms
9     FB Intercept         e  0.559687295  1.083121e-01 1.402673e+00   brms
10    FB Intercept        eh  1.791385030  7.392933e-01 4.180537e+00   brms
11    FB Intercept    r_tref  0.305160435  2.244636e-01 6.919619e-01   brms
12    FB Intercept        th 12.662596813  5.030311e+00 2.222652e+01   brms
13    FM Intercept         e  0.431295139 -1.453300e-01 1.333040e+00   brms
14    FM Intercept        eh  1.744530292  8.002445e-01 3.706611e+00   brms
15    FM Intercept    r_tref  0.292560194  2.058322e-01 6.685834e-01   brms
16    FM Intercept        th 12.043311409  4.883678e+00 2.118363e+01   brms
17    HO Intercept         e  0.709767004  1.675185e-01 1.623796e+00   brms
18    HO Intercept        eh  1.419402553  5.444447e-02 3.623925e+00   brms
19    HO Intercept    r_tref  0.260537905  1.585935e-01 5.697207e-01   brms
20    HO Intercept        th 11.430375290  4.122478e+00 2.131629e+01   brms
21    JM Intercept         e  0.524660276  5.999762e-02 1.379523e+00   brms
22    JM Intercept        eh  1.809604223  7.951560e-01 3.947093e+00   brms
23    JM Intercept    r_tref  0.297065772  2.147322e-01 6.788334e-01   brms
24    JM Intercept        th 12.325227102  4.997795e+00 2.139359e+01   brms
25    MU Intercept         e  0.499049714 -1.028082e-01 1.364656e+00   brms
26    MU Intercept        eh  1.684952381  6.300327e-01 3.962926e+00   brms
27    MU Intercept    r_tref  0.284802633  1.932996e-01 6.383678e-01   brms
28    MU Intercept        th 11.835542362  4.739944e+00 2.155500e+01   brms
29    RA Intercept         e  0.558131193  3.032568e-02 1.394866e+00   brms
30    RA Intercept        eh  1.650742366  5.117361e-01 3.949369e+00   brms
31    RA Intercept    r_tref  0.283150448  1.902202e-01 6.331985e-01   brms
32    RA Intercept        th 11.953127784  4.649346e+00 2.154704e+01   brms
33 SI_EK Intercept         e  0.515749586 -6.130383e-02 1.396251e+00   brms
34 SI_EK Intercept        eh  1.584378420  3.848152e-01 3.557932e+00   brms
35 SI_EK Intercept    r_tref  0.268151899  1.711672e-01 6.101190e-01   brms
36 SI_EK Intercept        th 11.291658804  4.389893e+00 2.093557e+01   brms
37 SI_HA Intercept         e  0.386536978 -2.911310e-01 1.333947e+00   brms
38 SI_HA Intercept        eh  1.629646899  7.275626e-01 3.175476e+00   brms
39 SI_HA Intercept    r_tref  0.281519344  1.873217e-01 6.468407e-01   brms
40 SI_HA Intercept        th 11.569018302  4.675350e+00 2.138970e+01   brms
41    JM      <NA>    r_tref  0.351191766 -7.483317e-01 1.450715e+00    nls
42    JM      <NA>         e  0.996605137 -2.595589e+00 4.588800e+00    nls
43    JM      <NA>        eh  2.386956961 -3.466243e-01 5.120538e+00    nls
44    JM      <NA>        th  9.696491298 -1.372360e+01 3.311658e+01    nls
45    FM      <NA>    r_tref  0.259069655 -2.973118e+00 3.491257e+00    nls
46    FM      <NA>         e  0.066736193 -8.449288e+00 8.582761e+00    nls
47    FM      <NA>        eh  1.168188267 -1.669096e+01 1.902734e+01    nls
48    FM      <NA>        th 16.854069553 -2.585282e+02 2.922363e+02    nls
49 SI_HA      <NA>    r_tref  0.351191766 -4.408465e+01 4.478703e+01    nls
50 SI_HA      <NA>         e  0.000000000 -5.795361e+01 5.795361e+01    nls
51 SI_HA      <NA>        eh  0.718116309 -1.020051e+01 1.163674e+01    nls
52 SI_HA      <NA>        th 10.516245975 -2.839112e+03 2.860144e+03    nls
53 SI_EK      <NA>    r_tref  0.192145129  1.245700e-01 2.597203e-01    nls
54 SI_EK      <NA>         e  0.869971502 -9.913036e-01 2.731247e+00    nls
55 SI_EK      <NA>        eh 13.511212154 -8.731800e+00 3.575422e+01    nls
56 SI_EK      <NA>        th  9.573238565  8.670834e+00 1.047564e+01    nls
57    FB      <NA>    r_tref  0.351191766 -1.207706e+02 1.214730e+02    nls
58    FB      <NA>         e  0.441722689 -1.144877e+02 1.153711e+02    nls
59    FB      <NA>        eh  0.669350807 -1.083993e+02 1.097380e+02    nls
60    FB      <NA>        th 14.912418904 -9.671099e+03 9.700923e+03    nls
61    BT      <NA>    r_tref  0.351191766 -1.090479e+00 1.792863e+00    nls
62    BT      <NA>         e  1.201395786 -4.225225e+00 6.628016e+00    nls
63    BT      <NA>        eh  1.938535110 -3.142411e-02 3.908494e+00    nls
64    BT      <NA>        th  9.568401144 -2.807795e+01 4.721475e+01    nls
65    MU      <NA>    r_tref  0.351191766 -3.367455e+04 3.367525e+04    nls
66    MU      <NA>         e  0.000000000 -9.691759e+03 9.691759e+03    nls
67    MU      <NA>        eh  0.206487598 -2.445567e+03 2.445980e+03    nls
68    MU      <NA>        th 14.144568959 -7.193905e+06 7.193934e+06    nls
69    HO      <NA>    r_tref  0.342950664 -3.632470e+02 3.639329e+02    nls
70    HO      <NA>         e  0.769380178 -7.131824e+02 7.147211e+02    nls
71    HO      <NA>        eh  0.000110422 -1.428046e+03 1.428046e+03    nls
72    HO      <NA>        th 16.854069553 -1.910990e+08 1.910990e+08    nls
73    BS      <NA>    r_tref  0.184145825  3.776306e-02 3.305286e-01    nls
74    BS      <NA>         e  1.340771288 -2.661068e+00 5.342611e+00    nls
75    BS      <NA>        eh  6.667623916 -1.819183e+01 3.152707e+01    nls
76    BS      <NA>        th 11.378692235  4.801183e+00 1.795620e+01    nls
77    RA      <NA>    r_tref  0.351191766 -8.840480e+03 8.841182e+03    nls
78    RA      <NA>         e  0.425716141 -3.266927e+03 3.267779e+03    nls
79    RA      <NA>        eh  0.271786041 -1.780731e+03 1.781275e+03    nls
80    RA      <NA>        th 16.854069553 -1.557304e+06 1.557337e+06    nls

Can we fit the brms model to areas separately and check the shrinkage (in coefficients and predictions) from using a hierarchical model?

dlist <- list()
dpred <- list()

for(i in unique(d$area)) {
  
  dd <- filter(dat, area == i)
  
  fit_area_spec <- brm(
    bf(rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))),
       rtref + e + eh + th ~ 1,
       nl = TRUE),
    data = dd,
    family = student,
    cores = 4,
      chains = 4,
      iter = 5000,
      #sample_prior = "yes",
    #save_pars = save_pars(all = TRUE),
    seed = 9,
    prior = prior,
    control = list(adapt_delta = 0.99, max_treedepth = 12)
)

  # Get coefficients
  brms_pars_area_spec <- fit_area_spec |> 
  gather_draws(b_rtref_Intercept,
               b_e_Intercept,
               b_eh_Intercept,
               b_th_Intercept) |> 
  rename(parameter = .variable) |> 
  mutate(parameter = str_remove(parameter, "_Intercept"),
         parameter = str_remove(parameter, "b_"),
         parameter = ifelse(parameter == "rtref", "r_tref", parameter)) |> 
  summarise(Estimate = median(.value),
            lwr = quantile(.value, probs = 0.025),
            upr = quantile(.value, probs = 0.975)) |> 
  ungroup() |> 
  mutate(source = "brms area spec.", 
         area = i)
  
  dlist[[i]] <- brms_pars_area_spec
  
  # Make predictions
  preds_brms <- dd |>
    data_grid(temp = seq_range(temp, n = 100)) |> 
    mutate(bk = 8.62e-05,
           tref = 8 + 273.15) |> 
    add_epred_draws(fit_area_spec, re_formula = NA) |> 
    ungroup() |>
    mutate(area = i)
  
  dpred[[i]] <- preds_brms

}

filter: removed 284 rows (93%), 22 rows remaining

Compiling Stan program...

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

Start sampling

rename: renamed one variable (parameter)

mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)

summarise: now 4 rows and 4 columns, ungrouped

ungroup: no grouping variables

mutate: new variable 'source' (character) with one unique value and 0% NA

        new variable 'area' (character) with one unique value and 0% NA

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

ungroup: no grouping variables

mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 289 rows (94%), 17 rows remaining

Compiling Stan program...

recompiling to avoid crashing R session

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

Start sampling

rename: renamed one variable (parameter)

mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)

summarise: now 4 rows and 4 columns, ungrouped

ungroup: no grouping variables

mutate: new variable 'source' (character) with one unique value and 0% NA

        new variable 'area' (character) with one unique value and 0% NA

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

ungroup: no grouping variables

mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 259 rows (85%), 47 rows remaining

Compiling Stan program...

recompiling to avoid crashing R session

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

Start sampling

rename: renamed one variable (parameter)

mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)

summarise: now 4 rows and 4 columns, ungrouped

ungroup: no grouping variables

mutate: new variable 'source' (character) with one unique value and 0% NA

        new variable 'area' (character) with one unique value and 0% NA

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

ungroup: no grouping variables

mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 269 rows (88%), 37 rows remaining

Compiling Stan program...

Start sampling

rename: renamed one variable (parameter)

mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)

summarise: now 4 rows and 4 columns, ungrouped

ungroup: no grouping variables

mutate: new variable 'source' (character) with one unique value and 0% NA

        new variable 'area' (character) with one unique value and 0% NA

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

ungroup: no grouping variables

mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 277 rows (91%), 29 rows remaining

Compiling Stan program...

Start sampling

rename: renamed one variable (parameter)

mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)

summarise: now 4 rows and 4 columns, ungrouped

ungroup: no grouping variables

mutate: new variable 'source' (character) with one unique value and 0% NA

        new variable 'area' (character) with one unique value and 0% NA

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

ungroup: no grouping variables

mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 246 rows (80%), 60 rows remaining

Compiling Stan program...

Start sampling

rename: renamed one variable (parameter)

mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)

summarise: now 4 rows and 4 columns, ungrouped

ungroup: no grouping variables

mutate: new variable 'source' (character) with one unique value and 0% NA

        new variable 'area' (character) with one unique value and 0% NA

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

ungroup: no grouping variables

mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 288 rows (94%), 18 rows remaining

Compiling Stan program...

Start sampling

rename: renamed one variable (parameter)

mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)

summarise: now 4 rows and 4 columns, ungrouped

ungroup: no grouping variables

mutate: new variable 'source' (character) with one unique value and 0% NA

        new variable 'area' (character) with one unique value and 0% NA

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

ungroup: no grouping variables

mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 292 rows (95%), 14 rows remaining

Compiling Stan program...

Start sampling

rename: renamed one variable (parameter)

mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)

summarise: now 4 rows and 4 columns, ungrouped

ungroup: no grouping variables

mutate: new variable 'source' (character) with one unique value and 0% NA

        new variable 'area' (character) with one unique value and 0% NA

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

ungroup: no grouping variables

mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 268 rows (88%), 38 rows remaining

Compiling Stan program...

Start sampling

rename: renamed one variable (parameter)

mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)

summarise: now 4 rows and 4 columns, ungrouped

ungroup: no grouping variables

mutate: new variable 'source' (character) with one unique value and 0% NA

        new variable 'area' (character) with one unique value and 0% NA

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

ungroup: no grouping variables

mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 282 rows (92%), 24 rows remaining

Compiling Stan program...

Start sampling

rename: renamed one variable (parameter)

mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)

summarise: now 4 rows and 4 columns, ungrouped

ungroup: no grouping variables

mutate: new variable 'source' (character) with one unique value and 0% NA

        new variable 'area' (character) with one unique value and 0% NA

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

ungroup: no grouping variables

mutate: new variable 'area' (character) with one unique value and 0% NA

# Save aic as data frames
area_spec <- bind_rows(dlist)
area_spec_preds <- bind_rows(dpred)

# Plot prediction comparison
ggplot(area_pred_brms, aes(temp, y = .epred, color = factor(area, order$area),
                           fill = factor(area, order$area))) +
  stat_lineribbon(.width = c(0.95), alpha = 0.1, color = NA) +
  stat_lineribbon(.width = c(0), alpha = 0.6) + 
  # add non-hierarchical model predictions
  stat_lineribbon(data = area_spec_preds, .width = 0, alpha = 0.6,
                  aes(temp, y = .epred, color = factor(area, order$area), linetype = "non-hierarchial")) + 
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 0.4, alpha = 0.4) +
  scale_color_manual(values = colors, name = "Area") +
  scale_fill_manual(values = colors, name = "Area") +
  scale_linetype_manual(values = 2) +
  facet_wrap(~factor(area, rev(order$area))) +
  guides(fill = "none",
         color = guide_legend(nrow = 1, reverse = TRUE, title.position = "top", title.hjust = 0.5,
                              override.aes = list(size = 1))) +
  theme(axis.title.y = ggtext::element_markdown(),
        legend.position = "bottom",
        legend.direction = "horizontal") +
  labs(x = "Temperature (°C)",
       y = "von Bertalanffy growth coefficient (*k*)",
       linetype = "")

# Now zoom in on RA and add the nls fit
ggplot(area_pred_brms |> filter(area == "RA"), aes(temp, y = .epred, color = factor(area, order$area),
                           fill = factor(area, order$area))) +
  stat_lineribbon(.width = c(0.95), alpha = 0.1, color = NA) +
  stat_lineribbon(.width = c(0), alpha = 0.6) + 
  # add non-hierarchical model predictions
  stat_lineribbon(data = area_spec_preds |> filter(area == "RA"), .width = 0, alpha = 0.6,
                  aes(temp, y = .epred, color = factor(area, order$area), linetype = "non-hierarchial")) + 
  # add nls
  geom_line(data = preds |> filter(area == "RA"), aes(temp, .fitted, linetype = "nls multstart"), color = "gray10") +
  geom_point(data = dat |> filter(area == "RA"), aes(temp, rate, color = factor(area, order$area)), size = 0.4, alpha = 0.4) +
  scale_color_manual(values = colors, name = "Area") +
  scale_fill_manual(values = colors, name = "Area") +
  scale_linetype_manual(values = c(2, 3)) +
  facet_wrap(~factor(area, rev(order$area))) +
  guides(fill = "none",
         color = guide_legend(nrow = 1, reverse = TRUE, title.position = "top", title.hjust = 0.5,
                              override.aes = list(size = 1))) +
  theme(axis.title.y = ggtext::element_markdown(),
        legend.position = "bottom",
        legend.direction = "horizontal") +
  labs(x = "Temperature (°C)",
       y = "von Bertalanffy growth coefficient (*k*)",
       linetype = "")

filter (grouped): removed 9,000,000 rows (90%), 1,000,000 rows remaining

filter: removed 9,000,000 rows (90%), 1,000,000 rows remaining

filter: removed 900 rows (90%), 100 rows remaining

filter: removed 292 rows (95%), 14 rows remaining

# Plot parameter comparison
brms_comp <- bind_rows(brms_pars, area_spec)

brms_comp |> 
  ggplot(aes(area, Estimate, color = source)) + 
  geom_point(position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0,
                position = position_dodge(width = 0.5)) + 
  coord_cartesian(ylim = c(-3, 20)) + 
  facet_wrap(~parameter, scales = "free")

Extra analysis

knitr::knit_exit()

Temperature in growing season instead of all year average

Fit Sharpe-Schoolfield model to K

By area

model <- 'sharpeschoolhigh_1981'

# Get starting values on full dataset for Sharpe-Schoolfield model
dat2 <- VBG |>
  select(k_median, temp_gs, area) |>
  rename(rate = k_median,
         temp = temp_gs) # growing season! not annual average!

lower <- get_lower_lims(dat2$temp, dat2$rate, model_name = model)
upper <- get_upper_lims(dat2$temp, dat2$rate, model_name = model)
start <- get_start_vals(dat2$temp, dat2$rate, model_name = model)
  
# Sharpe-Schoolfield model fit to data for each area
preds2 <- NULL
pred2 <- NULL

for(a in unique(dat2$area)) {
  
  # Get data
  dd <- dat2 |> filter(area == a)
  
  # Fit model
  fit2 <- nls_multstart(
    rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 8),
    data = dd,
    iter = c(3, 3, 3, 3),
    start_lower = start*0.5,
    start_upper = start*2,
    lower = lower,
    upper = upper,
    supp_errors = 'Y'
    )
  
  # Make predictions on new data
  new_data2 <- data.frame(temp = seq(min(dd$temp), max(dd$temp), length.out = 100))
  
  pred2 <- augment(fit2, newdata = new_data2) |> mutate(area = a)
  
  # Add to general data frame
  preds2 <- data.frame(rbind(preds2, pred2))
  
}

All areas pooled

# One sharpe schoolfield fitted to all data
fit_all2 <- nls_multstart(
    rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 8),
    data = dat2,
    iter = c(3, 3, 3, 3),
    start_lower = start*0.5,
    start_upper = start*2,
    lower = lower,
    upper = upper,
    supp_errors = 'Y'
    )

# Make predictions on new data
new_data_all2 <- data.frame(temp = seq(min(dat2$temp), max(dat2$temp), length.out = 100))

pred_all2 <- augment(fit_all2, newdata = new_data_all2) |> 
  mutate(area = "all")

# Add t_opt
kb <- 8.62e-05
data.frame(par = names(coef(fit_all2)), est = coef(fit_all2)) |> 
  pivot_wider(names_from = par, values_from = est) |> 
  summarise(t_opt = (eh*th) / (eh + kb*th*log((eh/e)-1)))

Plot Sharpe Schoolfield fits

# Plot growth coefficients by year and area against mean SST
preds2 |>
  ggplot(aes(temp, .fitted, color = factor(area, order$area))) + 
  geom_point(data = dat2, aes(temp, rate, color = factor(area, order$area)), size = 0.6) +
  geom_line(linewidth = 1) +
  geom_line(data = pred_all2, aes(temp, .fitted), linewidth = 1, inherit.aes = FALSE, linetype = 2) +
  scale_color_manual(values = colors, name = "Area") +
  guides(color = guide_legend(override.aes = list(size = 1))) +
  scale_x_continuous(breaks = seq(-5, 20, 1)) +
  labs(x = "Temperature (°C)",
       y = "von Bertalanffy growth coefficient (*k*)") +
  theme(axis.title.y = ggtext::element_markdown())

Linf instead of k

By area

model <- 'sharpeschoolhigh_1981'

# Get starting values on full dataset for Sharpe-Schoolfield model
dat3 <- VBG |>
  select(linf_median, temp, area) |>
  rename(rate = linf_median)

lower <- get_lower_lims(dat3$temp, dat3$rate, model_name = model)
upper <- get_upper_lims(dat3$temp, dat3$rate, model_name = model)
start <- get_start_vals(dat3$temp, dat3$rate, model_name = model)
  
# Sharpe-Schoolfield model fit to data for each area
preds3 <- NULL
pred3 <- NULL

for(a in unique(dat3$area)) {
  
  # Get data
  dd <- dat3 |> filter(area == a)
  
  # Fit model
  fit3 <- nls_multstart(
    rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 8),
    data = dd,
    iter = c(3, 3, 3, 3),
    start_lower = start*0.5,
    start_upper = start*2,
    lower = lower,
    upper = upper,
    supp_errors = 'Y'
    )
  
  # Make predictions on new data
  new_data3 <- data.frame(temp = seq(min(dd$temp), max(dd$temp), length.out = 100))
  
  pred3 <- augment(fit3, newdata = new_data3) |> mutate(area = a)
  
  # Add to general data frame
  preds3 <- data.frame(rbind(preds3, pred3))
  
}

All areas pooled

# One sharpe schoolfield fitted to all data
fit_all3 <- nls_multstart(
    rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 8),
    data = dat3,
    iter = c(3, 3, 3, 3),
    start_lower = start*0.5,
    start_upper = start*2,
    lower = lower,
    upper = upper,
    supp_errors = 'Y'
    )

# Make predictions on new data
new_data_all3 <- data.frame(temp = seq(min(dat3$temp), max(dat3$temp), length.out = 100))

pred_all3 <- augment(fit_all3, newdata = new_data_all3) |> 
  mutate(area = "all")

Plot Sharpe Schoolfield fits

# Plot growth coefficients by year and area against mean SST
preds3 |>
  ggplot(aes(temp, .fitted, color = factor(area, order$area))) + 
  geom_point(data = dat3, aes(temp, rate, color = factor(area, order$area)), size = 0.6) +
  geom_line(linewidth = 1) +
  geom_line(data = pred_all3, aes(temp, .fitted), linewidth = 1, inherit.aes = FALSE, linetype = 2) +
  scale_color_manual(values = colors, name = "Area") +
  guides(color = guide_legend(override.aes = list(size = 1))) +
  scale_x_continuous(breaks = seq(-5, 20, 1)) +
  labs(x = "Temperature (°C)",
       y = expression(paste(italic(L[infinity]), " [cm]"))) +
  theme(axis.title.y = ggtext::element_markdown())

ggplot(VBG, aes(linf_median, k_median)) + 
  geom_point()